Variation in vocal behavior under varying stress levels
Budgie INBRE stress and contact calls
Marcelo Araya-Salas, PhD & Timothy Wright, PhD
2023-02-14
Source code and data found at https://github.com/maRce10/budgie_inbre_stress_vocal_output
Load packages
## add 'developer' to packages to be installed from github
x <- c("remotes", "lubridate", "readxl", "pbapply", "viridis", "ggplot2",
"kableExtra", "knitr", "formatR", "MASS", "sp", "GGally", "brms",
"lme4", "dplyr", "purrr", "forcats", "tidyr", "modelr", "tidybayes",
"cowplot", "ggrepel", "posterior", "ggridges", "maRce10/PhenotypeSpace")
source("~/Dropbox/R_package_testing/sketchy/R/load_packages.R")
load_packages(x)
source("~/Dropbox/R_package_testing/brmsish/R/html_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")
source("~/Dropbox/R_package_testing/brmsish/R/read_summary.R")Functions and global parameters
opts_knit$set(root.dir = "..")
# set evaluation false
opts_chunk$set(fig.width = 10, fig.height = 6, warning = FALSE, message = FALSE,
tidy = TRUE)
read_excel_df <- function(...) data.frame(read_excel(...))
# for reading months in english format
sl <- Sys.setlocale(locale = "en_US.UTF-8")
standard_error <- function(x) sd(x)/sqrt(length(x))
cols <- viridis(10, alpha = 0.7)
col_pointrange <- cols[7]# read ext sel tab calls
sels <- read.csv("./data/processed/tailored_budgie_calls_sel_tab.csv")
# keep only spectrographic parameters
sels <- sels[, c("sound.files", "selec", "duration", "meanfreq", "sd",
"freq.median", "freq.IQR", "time.IQR", "skew", "kurt", "sp.ent",
"time.ent", "entropy", "meandom", "mindom", "maxdom", "dfrange",
"modindx", "meanpeakf")]
sels$ID <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][1])
sels$month <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][2])
sels$day <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][3])
sels$year <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][4])
sels$date <- paste(sels$day, substr(sels$month, 0, 3), sels$year,
sep = "-")
sels$date <- as.Date(sels$date, format = "%d-%b-%Y")
# acoustic measurements
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")
indiv_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_first_week_by_individual.RDS")
indiv_ovlp$treatment <- factor(indiv_ovlp$treatment, levels = c("Control",
"Medium Stress", "High Stress"))
group_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")
# group_ovlp <- do.call(rbind, lapply(unique(group_ovlp$ID),
# function(x) { X <- group_ovlp[group_ovlp$ID == x, ]
# X$overlap.to.group <- X$overlap.to.group -
# X$overlap.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] X$distance.to.group <-
# X$distance.to.group - X$distance.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] return(X) }))
group_ovlp$treatment <- factor(group_ovlp$treatment, levels = c("Control",
"Medium Stress", "High Stress"))Counts per individual
df <- as.data.frame(table(sels$ID))
names(df) <- c("ID", "Sample_size")
df <- df[order(df$Sample_size, decreasing = FALSE), ]
kb <- kable(df, row.names = FALSE)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover",
"condensed", "responsive"))
print(kb)| ID | Sample_size |
|---|---|
| 132YGMM | 6 |
| 125YGMM | 12 |
| 160YGHM | 16 |
| 310YGHM | 21 |
| 393YGHM | 25 |
| 303YGHM | 37 |
| 398WBHM | 41 |
| 365YGHM | 44 |
| 399YGLM | 46 |
| 300YGHM | 47 |
| 270YGHM | 64 |
| 407YGHM | 97 |
| 386WBHM | 100 |
| 377WWLM | 108 |
| 367WWNM | 119 |
| 371YYLM | 149 |
| 397YGHM | 175 |
| 378YYLM | 195 |
| 404WBHM | 196 |
| 258YGHM | 223 |
| 389WWLM | 230 |
| 262YGHM | 306 |
| 400YGHM | 360 |
| 388YGLM | 377 |
| 327YYLM | 404 |
| 396YBHM | 455 |
| 403WBHM | 566 |
| 356WBHM | 761 |
| 361YGHM | 767 |
| 402YGHM | 849 |
| 395WBHM | 854 |
| 360YGHM | 900 |
| 390WBHM | 935 |
| 405YBHM | 975 |
| 385YBHM | 1256 |
| 394YBHM | 1693 |
Add metadata
metadat <- read_excel_df("./data/raw/INBREStress_MasterDataSheet_14Nov19.xlsx")
# head(metadat)
sels$ID[sels$ID == "125YGMM"] <- "125YGHM"
sels$ID[sels$ID == "394YBHM"] <- "394WBHM"
# setdiff(sels$ID, metadat$Bird.ID) setdiff(metadat$Bird.ID,
# sels$ID)
sels$treatment <- sapply(1:nrow(sels), function(x) {
metadat$Treatment[metadat$Bird.ID == sels$ID[x]][1]
})
sels$treatment.room <- sapply(1:nrow(sels), function(x) {
metadat$Treatment.Room[metadat$Bird.ID == sels$ID[x]][1]
})
sels$round <- sapply(1:nrow(sels), function(x) {
metadat$Round[metadat$Bird.ID == sels$ID[x]][1]
})
sels$source.room <- sapply(1:nrow(sels), function(x) {
metadat$Source.Room[metadat$Bird.ID == sels$ID[x]][1]
})
sels$record.group <- sapply(1:nrow(sels), function(x) {
metadat$Record.Group[metadat$Bird.ID == sels$ID[x]][1]
})
# add week
out <- lapply(unique(sels$round), function(x) {
Y <- sels[sels$round == x, ]
min_date <- min(Y$date)
week_limits <- min_date + seq(0, 100, by = 7)
Y$week <- NA
for (i in 2:length(week_limits)) Y$week[Y$date >= week_limits[i -
1] & Y$date < week_limits[i]] <- i - 1
return(Y)
})
sels <- do.call(rbind, out)
sels$cort.baseline <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$cort.stress <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$stress.response <- sels$cort.stress #- sels$cort.baseline
sels$weight <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$breath.count <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
# remove week 5
sels <- sels[sels$week != 5, ]agg_dat <- aggregate(selec ~ ID + week, data = sels, length)
# compare to week 1 agg_dat$call.count <-
# sapply(1:nrow(agg_dat), function(x) { baseline <-
# agg_dat$selec[agg_dat$week == 1 & agg_dat$ID == agg_dat$ID[x]]
# if (length(baseline) > 0) change <- agg_dat$selec[x] -
# baseline else change <- agg_dat$selec[x] return(change) } )
# without comparing to week 1
agg_dat$call.count <- sapply(1:nrow(agg_dat), function(x) agg_dat$selec[x])
agg_dat$selec <- NULL
agg_dat$baseline.CORT <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$cort.baseline[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$cort.baseline[sels$week == agg_dat$week[x] & sels$ID ==
agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$stress.response <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$stress.response[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$stress.response[sels$week == agg_dat$week[x] &
sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$stress.CORT <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$cort.stress[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$cort.stress[sels$week == agg_dat$week[x] & sels$ID ==
agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$weight <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$weight[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$weight[sels$week == agg_dat$week[x] & sels$ID ==
agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$breath.rate <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$breath.count[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$breath.count[sels$week == agg_dat$week[x] & sels$ID ==
agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$acoustic.diversity <- sapply(1:nrow(agg_dat), function(x) {
area <- areas_by_week$raref.area[areas_by_week$ID == agg_dat$ID[x] &
areas_by_week$week == agg_dat$week[x]]
if (length(area) < 1)
area <- NA
return(area)
})
agg_dat$acoustic.distance <- sapply(1:nrow(agg_dat), function(x) {
distance <- indiv_ovlp$distance.to.first.week[indiv_ovlp$ID ==
agg_dat$ID[x] & indiv_ovlp$week == agg_dat$week[x]]
if (length(distance) < 1)
distance <- NA
return(distance)
})
agg_dat$acustic.plasticity <- sapply(1:nrow(agg_dat), function(x) {
overlap <- indiv_ovlp$overlap.to.first.week[indiv_ovlp$ID == agg_dat$ID[x] &
indiv_ovlp$week == agg_dat$week[x]]
plasticity <- 1 - overlap
if (length(plasticity) < 1)
plasticity <- NA
return(plasticity)
})
agg_dat$acoustic.convergence <- sapply(1:nrow(agg_dat), function(x) {
overlap <- group_ovlp$overlap.to.group[group_ovlp$ID == agg_dat$ID[x] &
group_ovlp$week == agg_dat$week[x]]
if (length(overlap) < 1)
overlap <- NA
return(overlap)
})
agg_dat$treatment <- sapply(1:nrow(agg_dat), function(x) unique(sels$treatment[sels$ID ==
agg_dat$ID[x]]))
agg_dat$round <- sapply(1:nrow(agg_dat), function(x) unique(sels$round[sels$ID ==
agg_dat$ID[x]]))Physiological parameters
Barplot and effect sizes graph
physio_models <- readRDS("./data/processed/physiological_response_models.RDS")
breath.count <- stack(metadat[, c("D3.Breath.Count", "D7.Breath.Count",
"D14.Breath.Count", "D21.Breath.Count", "D28.Breath.Count")])
weight <- stack(metadat[, c("D3.Bird.Weight..g.", "D7.Bird.Weight..g.",
"D14.Bird.Weight..g.", "D21.Bird.Weight..g.", "D28.Bird.Weight..g.")])
cort.stress <- stack(metadat[, c("D3.CORT.Stress", "D7.CORT.Stress",
"D14.CORT.Stress", "D21.CORT.Stress", "D28.CORT.Stress")])
cort.baseline <- stack(metadat[, c("D3.CORT.Baseline", "D7.CORT.Baseline",
"D14.CORT.Baseline", "D21.CORT.Baseline", "D28.CORT.Baseline")])
stress <- data.frame(metadat[, c("Bird.ID", "Treatment", "Round",
"Treatment.Room")], week = breath.count$ind, breath.count = breath.count$values,
weight = weight$values, cort.stress = cort.stress$values, cort.baseline = cort.baseline$values,
stress.response = cort.stress$values - cort.baseline$values)
# head(stress)
stress$week <- factor(sapply(strsplit(as.character(stress$week), "\\."),
"[[", 1), levels = c("D3", "D7", "D14", "D21", "D28"))
stress$day <- as.numeric(gsub("D", "", as.character(stress$week)))
stress$round <- paste("Round", stress$Round)
stress$treatment <- factor(stress$Treatment, levels = c("Control",
"Medium Stress", "High Stress"))
# remove 5th week
stress <- stress[stress$week != "D28", ]
stress_l <- lapply(stress$Bird.ID, function(x) {
X <- stress[stress$Bird.ID == x, ]
X$breath.count <- X$breath.count - X$breath.count[X$week == "D3"]
X$weight <- X$weight - X$weight[X$week == "D3"]
X$cort.stress <- X$cort.stress - X$cort.stress[X$week == "D3"]
X$cort.baseline <- X$cort.baseline - X$cort.baseline[X$week ==
"D3"]
X$stress.response <- X$stress.response #- X$stress.response[X$week == 'D3']
return(X)
})
stress <- do.call(rbind, stress_l)
agg_stress <- aggregate(cbind(breath.count, weight, stress.response,
cort.baseline) ~ week + treatment, stress, mean)
agg_stress_se <- aggregate(cbind(breath.count, weight, stress.response,
cort.baseline) ~ week + treatment, stress, standard_error)
names(agg_stress_se) <- paste(names(agg_stress_se), ".se", sep = "")
agg_stress <- cbind(agg_stress, agg_stress_se[, 3:6])
agg_stress$Week <- 1:4
bs <- 10
gg_breath.count <- ggplot(data = agg_stress, aes(x = Week, y = breath.count,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = breath.count -
breath.count.se, ymax = breath.count + breath.count.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Mean change in breath\nrate (breaths/min)",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_weight <- ggplot(data = agg_stress, aes(x = Week, y = weight, fill = treatment)) +
geom_bar(stat = "identity") + geom_errorbar(aes(ymin = weight -
weight.se, ymax = weight + weight.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") +
labs(y = "Mean change in \nweight (grams)", x = "Week") + theme_classic(base_size = bs) +
theme(legend.position = "none")
gg_cort.baseline <- ggplot(data = agg_stress, aes(x = Week, y = cort.baseline,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = cort.baseline -
cort.baseline.se, ymax = cort.baseline + cort.baseline.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Mean change in\nbaseline CORT (ng/mL)",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_stress.response <- ggplot(data = agg_stress, aes(x = Week, y = stress.response,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = stress.response -
stress.response.se, ymax = stress.response + stress.response.se),
width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Mean change in stress response \nCORT (ng/mL)",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_coeffs_physio <- lapply(physio_models, function(x) {
vars <- grep("b_", posterior::variables(x), value = TRUE)
draws <- posterior::as_draws_array(x, variable = vars)
coef_table <- draw_summary(draws, variables = vars, probs = c(0.025,
0.975), robust = TRUE)
coef_table$predictor <- rownames(coef_table)
coef_table$predictor <- gsub("b_treatment|b_", "", coef_table$predictor)
coef_table$predictor <- gsub("Stress", " stress", coef_table$predictor)
coef_table$predictor <- gsub("week", "Week", coef_table$predictor)
coef_table <- coef_table[coef_table$predictor != "Intercept",
]
gg_coef <- ggplot2::ggplot(data = coef_table, aes(x = Estimate,
y = predictor)) + geom_vline(xintercept = 0, lty = 2) + ggplot2::geom_point(size = 4,
col = col_pointrange) + ggplot2::geom_errorbar(ggplot2::aes(xmin = `l-95% CI`,
xmax = `u-95% CI`), width = 0, col = col_pointrange) + ggplot2::theme_classic(base_size = bs) +
ggplot2::theme(axis.ticks.length = ggplot2::unit(0, "pt"),
plot.margin = ggplot2::margin(0, 0, 0, 0, "pt"), legend.position = "none",
strip.background = ggplot2::element_blank(), strip.text = ggplot2::element_blank()) +
ggplot2::labs(x = "Effect size", y = "")
return(gg_coef)
})
cowplot::plot_grid(gg_breath.count, gg_weight, gg_cort.baseline, gg_stress.response,
gg_coeffs_physio[[grep("breath", names(gg_coeffs_physio))]] +
theme_classic(base_size = bs), gg_coeffs_physio[[grep("weight",
names(gg_coeffs_physio))]] + theme_classic(base_size = bs),
gg_coeffs_physio[[grep("baseline", names(gg_coeffs_physio))]] +
theme_classic(base_size = bs), gg_coeffs_physio[[grep("response",
names(gg_coeffs_physio))]] + theme_classic(base_size = bs),
nrow = 2, rel_heights = c(1.8, 1))
# try bs = 20 for saving plots
# cowplot::ggsave2(filename =
# './output/bar_graphs_and_estimates_physiology_70dpi.jpeg',
# width = 25, height = 9) cowplot::ggsave2(filename =
# './output/bar_graphs_and_estimates_physiology_300dpi.jpeg',
# dpi = 300, width = 25, height = 9)Stats
Models: Predicted physio measure ~ treatment + week (continuous) + IndRandom
Variables (Difference from Week 1): weight, BR, baseline CORT, Stress CORT, Stress Response
responses <- c("baseline.CORT", "stress.response", "stress.CORT",
"weight", "breath.rate")
predictors <- c("~ treatment + week + (1|ID) + (1|round)")
formulas <- expand.grid(responses = responses, predictors = predictors,
stringsAsFactors = FALSE)
vars_to_scale <- c(responses, "week")
# remove week 1
sub_agg_dat <- agg_dat[agg_dat$week != 1, ]
for (i in vars_to_scale) sub_agg_dat[, vars_to_scale] <- scale(sub_agg_dat[,
vars_to_scale])
physio_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- sub_agg_dat[!is.na(sub_agg_dat[names(sub_agg_dat) ==
formulas$responses[x]]), ]
sub_dat
mod <- brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 20000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
chains = 4, prior = c(prior(normal(0, 5), "b"), prior(normal(0,
10), "Intercept"), prior(student_t(3, 0, 10), "sd"), prior(student_t(3,
0, 10), "sigma")))
return(mod)
})
names(physio_models) <- paste(formulas$responses, formulas$predictors)
saveRDS(physio_models, "./data/processed/physiological_response_models.RDS")physio_models <- readRDS("./data/processed/physiological_response_models.RDS")
ggl <- list()
for (x in 1:length(physio_models)) {
ggl[[length(ggl) + 1]] <- html_summary(physio_models[[x]], gsub.pattern = "b_treatment|b_",
gsub.replacement = "", highlight = FALSE, remove.intercepts = TRUE,
model.name = names(physio_models)[x])
print(ggl[[length(ggl)]])
}baseline.CORT ~ treatment + week + (1|ID) + (1|round)
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | baseline.CORT ~ treatment + week + (1 | ID) + (1 | round) | 20000 | 4 | 1 | 10000 | 406 | 0 | 12756.08 | 17437.15 | 614553873 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | 0.760 | 0.108 | 1.420 | 1 | 13169.12 | 17881.30 |
| MediumStress | 0.180 | -0.506 | 0.868 | 1 | 12756.08 | 17437.15 |
| week | -0.001 | -0.142 | 0.141 | 1 | 27598.67 | 26244.08 |
stress.response ~ treatment + week + (1|ID) + (1|round)
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | stress.response ~ treatment + week + (1 | ID) + (1 | round) | 20000 | 4 | 1 | 10000 | 2129 | 0 | 1831.829 | 3145.311 | 1852853265 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | -0.204 | -0.889 | 0.465 | 1.003 | 1831.829 | 3145.311 |
| MediumStress | 0.097 | -0.667 | 0.839 | 1.002 | 2582.692 | 10967.477 |
| week | -0.152 | -0.284 | -0.020 | 1.003 | 6379.514 | 15589.581 |
stress.CORT ~ treatment + week + (1|ID) + (1|round)
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | stress.CORT ~ treatment + week + (1 | ID) + (1 | round) | 20000 | 4 | 1 | 10000 | 1015 | 0 | 7696.824 | 8195.241 | 672971911 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | -0.219 | -0.905 | 0.452 | 1.001 | 8677.295 | 16469.149 |
| MediumStress | 0.084 | -0.684 | 0.832 | 1 | 7696.824 | 8195.241 |
| week | -0.154 | -0.284 | -0.022 | 1 | 25064.631 | 25062.987 |
weight ~ treatment + week + (1|ID) + (1|round)
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | weight ~ treatment + week + (1 | ID) + (1 | round) | 20000 | 4 | 1 | 10000 | 1142 | 0 | 6459.84 | 11029.72 | 351999546 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | -0.319 | -0.971 | 0.328 | 1.001 | 10517.09 | 16418.60 |
| MediumStress | -0.120 | -0.850 | 0.598 | 1.002 | 6459.84 | 11029.72 |
| week | -0.345 | -0.478 | -0.211 | 1 | 14142.17 | 16891.28 |
breath.rate ~ treatment + week + (1|ID) + (1|round)
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | breath.rate ~ treatment + week + (1 | ID) + (1 | round) | 20000 | 4 | 1 | 10000 | 732 | 0 | 17760.6 | 14209.36 | 1444112821 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | 0.147 | -0.156 | 0.449 | 1 | 17760.60 | 14209.36 |
| MediumStress | 0.060 | -0.264 | 0.389 | 1 | 18526.45 | 24064.77 |
| week | -0.784 | -0.898 | -0.670 | 1 | 24048.27 | 26983.78 |
names(ggl) <- names(physio_models)
Takeaways
Breath rate decreases gradually with time across after the first week
Stress response is higher in “high stress” birds compared to first week
Acoustic space projection
t-SNE
scale_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median",
"freq.IQR", "time.IQR", "skew", "kurt", "sp.ent", "time.ent",
"entropy", "meandom", "mindom", "maxdom", "dfrange", "modindx",
"meanpeakf")])
tsne <- Rtsne(scale_param, dims = 2, perplexity = 30, verbose = FALSE,
max_iter = 5000)
saveRDS(tsne, "./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")tsne <- readRDS("./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")
Y <- as.data.frame(tsne$Y)
names(Y) <- c("TSNE1", "TSNE2")
sels <- data.frame(sels, Y)
sels$treatment <- factor(sels$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
ggplot(sels, aes(x = TSNE1, y = TSNE2, col = as.factor(treatment))) +
geom_point() + labs(color = "Treatment") + scale_color_viridis_d(alpha = 0.4) +
theme_classic(base_size = 25) + guides(colour = guide_legend(override.aes = list(size = 10)))
Behavioral parameters
Barplot and effect sizes pgrah
behav_models <- readRDS("./data/processed/behavioral_response_models.RDS")
agg_call.count <- aggregate(cbind(call.count, acoustic.convergence) ~
week + treatment, agg_dat, mean)
agg_behav <- aggregate(cbind(acoustic.diversity, acustic.plasticity) ~
week + treatment, agg_dat, mean)
agg_call.count_se <- aggregate(cbind(call.count, acoustic.convergence) ~
week + treatment, agg_dat, standard_error)
agg_behav_se <- aggregate(cbind(acoustic.diversity, acustic.plasticity) ~
week + treatment, agg_dat, standard_error)
agg_behav_se <- merge(agg_call.count_se, agg_behav_se, all = TRUE)
names(agg_behav_se) <- paste(names(agg_behav_se), ".se", sep = "")
agg_behav <- merge(agg_call.count, agg_behav, all = TRUE)
agg_behav <- cbind(agg_behav, agg_behav_se[, 3:6])
bs <- 10
agg_behav$treatment <- factor(agg_behav$treatment, levels = c("Control",
"Medium Stress", "High Stress"))
gg_call.count <- ggplot(data = agg_behav, aes(x = week, y = call.count,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = call.count -
call.count.se, ymax = call.count + call.count.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Vocal output", x = "Week") +
theme_classic(base_size = bs) + theme(legend.position = "none")
gg_acoustic.diversity <- ggplot(data = agg_behav, aes(x = week, y = acoustic.diversity,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acoustic.diversity -
acoustic.diversity.se, ymax = acoustic.diversity + acoustic.diversity.se),
width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Change in vocal diversity",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_acustic.plasticity <- ggplot(data = agg_behav, aes(x = week, y = acustic.plasticity,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acustic.plasticity -
acustic.plasticity.se, ymax = acustic.plasticity + acustic.plasticity.se),
width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Vocal plasticity",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_acoustic.convergence <- ggplot(data = agg_behav, aes(x = week,
y = acoustic.convergence, fill = treatment)) + geom_bar(stat = "identity") +
geom_errorbar(aes(ymin = acoustic.convergence - acoustic.convergence.se,
ymax = acoustic.convergence + acoustic.convergence.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Vocal convergence", x = "Week") +
theme_classic(base_size = bs) + theme(legend.position = "none")
gg_coeffs_behav <- lapply(behav_models, function(x) {
vars <- grep("b_", posterior::variables(x), value = TRUE)
draws <- posterior::as_draws_array(x, variable = vars)
coef_table <- draw_summary(draws, variables = vars, probs = c(0.025,
0.975), robust = TRUE)
coef_table$predictor <- rownames(coef_table)
coef_table$predictor <- gsub("b_treatment|b_", "", coef_table$predictor)
coef_table$predictor <- gsub("Stress", " stress", coef_table$predictor)
coef_table$predictor <- gsub("week", "Week", coef_table$predictor)
coef_table <- coef_table[coef_table$predictor != "Intercept",
]
gg_coef <- ggplot2::ggplot(data = coef_table, aes(x = Estimate,
y = predictor)) + geom_vline(xintercept = 0, lty = 2) + ggplot2::geom_point(size = 4,
col = col_pointrange) + ggplot2::geom_errorbar(ggplot2::aes(xmin = `l-95% CI`,
xmax = `u-95% CI`), width = 0, col = col_pointrange) + ggplot2::theme_classic(base_size = bs) +
ggplot2::theme(axis.ticks.length = ggplot2::unit(0, "pt"),
plot.margin = ggplot2::margin(0, 0, 0, 0, "pt"), legend.position = "none",
strip.background = ggplot2::element_blank(), strip.text = ggplot2::element_blank()) +
ggplot2::labs(x = "Effect size", y = "")
return(gg_coef)
})
cowplot::plot_grid(gg_call.count, gg_acoustic.diversity, gg_acustic.plasticity,
gg_acoustic.convergence, gg_coeffs_behav[[grep("count", names(gg_coeffs_behav))]] +
theme_classic(base_size = bs), gg_coeffs_behav[[grep("diversity",
names(gg_coeffs_behav))]] + theme_classic(base_size = bs),
gg_coeffs_behav[[grep("plasticity", names(gg_coeffs_behav))]] +
theme_classic(base_size = bs), gg_coeffs_behav[[grep("convergence",
names(gg_coeffs_behav))]] + theme_classic(base_size = bs),
nrow = 2, rel_heights = c(1.8, 1))
# try bs = 20 for saving plots
# cowplot::ggsave2(filename =
# './output/bar_graphs_and_estimates_behavior_70dpi.jpeg', width
# = 25, height = 9) # cowplot::ggsave2(filename =
# './output/bar_graphs_and_estimates_behavior_300dpi.jpeg', dpi
# = 300, width = 25, height = 9)Stats
Model: Predicted behavior ~ treatment + week (continuous) + IndRandom
Variables: # calls, Distance moved from self in first week, Overlap to original acoustic space, Match to group repertoire, Maybe overall size of acoustic space
Do as comparison to week one using rarefacted calls and kernel density
responses <- c("call.count", "acoustic.diversity", "acustic.plasticity",
"acoustic.convergence")
predictors <- c("~ treatment + week + (1|ID) + (1|round)")
formulas <- expand.grid(responses = responses, predictors = predictors,
stringsAsFactors = FALSE)
vars_to_scale <- c(responses, "week")
for (i in vars_to_scale) agg_dat[, vars_to_scale] <- scale(agg_dat[,
vars_to_scale])
behav_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
]
# remove week 1
if (!grepl("count|group", formulas$responses[x]))
sub_dat <- sub_dat[sub_dat$week != 1, ]
mod <- brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 50000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9,
max_treedepth = 15), chains = 4, prior = c(prior(normal(0,
5), "b"), prior(normal(0, 10), "Intercept"), prior(student_t(3,
0, 10), "sd"), prior(student_t(3, 0, 10), "sigma")))
return(mod)
})
names(behav_models) <- paste(formulas$responses, formulas$predictors)
saveRDS(behav_models, "./data/processed/behavioral_response_models.RDS")behav_models <- readRDS("./data/processed/behavioral_response_models.RDS")
ggl <- list()
for (x in 1:length(behav_models)) {
ggl[[length(ggl) + 1]] <- html_summary(behav_models[[x]], gsub.pattern = "b_treatment|b_",
gsub.replacement = "", highlight = FALSE, remove.intercepts = TRUE,
model.name = names(behav_models)[x])
print(ggl[[length(ggl)]])
}call.count ~ treatment + week + (1|ID) + (1|round)
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | call.count ~ treatment + week + (1 | ID) + (1 | round) | 50000 | 4 | 1 | 25000 | 1455 | 0 | 21916.29 | 19215.21 | 2086259884 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | 0.784 | 0.205 | 1.357 | 1 | 21916.29 | 19215.21 |
| MediumStress | 0.394 | -0.183 | 0.971 | 1 | 29122.91 | 29122.87 |
| week | -0.059 | -0.203 | 0.086 | 1 | 48490.61 | 55075.00 |
acoustic.diversity ~ treatment + week + (1|ID) + (1|round)
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | acoustic.diversity ~ treatment + week + (1 | ID) + (1 | round) | 50000 | 4 | 1 | 25000 | 2953 | 0 | 6138.296 | 3258.34 | 1000749341 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | -0.208 | -1.070 | 0.657 | 1.001 | 6138.296 | 3258.34 |
| MediumStress | -0.699 | -1.559 | 0.158 | 1.001 | 8701.168 | 12053.83 |
| week | 0.036 | -0.081 | 0.152 | 1 | 27120.995 | 32683.11 |
acustic.plasticity ~ treatment + week + (1|ID) + (1|round)
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | acustic.plasticity ~ treatment + week + (1 | ID) + (1 | round) | 50000 | 4 | 1 | 25000 | 7397 | 0 | 1025.75 | 483.942 | 803966074 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | -0.904 | -1.731 | -0.062 | 1.003 | 4966.382 | 30252.390 |
| MediumStress | -0.232 | -1.130 | 0.678 | 1.005 | 1025.750 | 483.942 |
| week | 0.171 | 0.034 | 0.306 | 1.002 | 2427.679 | 26782.955 |
acoustic.convergence ~ treatment + week + (1|ID) + (1|round)
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | acoustic.convergence ~ treatment + week + (1 | ID) + (1 | round) | 50000 | 4 | 1 | 25000 | 8224 | 0 | 150.4 | 57.174 | 1972081133 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | 0.335 | -0.307 | 1.003 | 1.006 | 1086.629 | 21463.246 |
| MediumStress | 0.372 | -0.288 | 1.045 | 1.003 | 12829.194 | 16680.338 |
| week | -0.228 | -0.406 | -0.058 | 1.017 | 150.400 | 57.174 |
names(ggl) <- names(behav_models)
Takeaways
Higher vocal output in “high stress” birds compared to control
Higher acoustic overlap to themselves in week 1 for “high stress” birds compared to control
Decrease in overlap to themselves in week 1 across time
Session information
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=pt_BR.UTF-8
## [7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] PhenotypeSpace_0.1.0 warbleR_1.1.26 NatureSounds_1.0.4
## [4] seewave_2.1.8 tuneR_1.3.3.1 ggridges_0.5.3
## [7] posterior_1.2.0 ggrepel_0.9.1 cowplot_1.1.1
## [10] tidybayes_3.0.2 modelr_0.1.8 tidyr_1.2.0
## [13] forcats_0.5.1 purrr_0.3.4 dplyr_1.0.8
## [16] lme4_1.1-28 Matrix_1.3-4 brms_2.16.3
## [19] Rcpp_1.0.8 GGally_2.1.2 sp_1.4-6
## [22] MASS_7.3-54 formatR_1.11 knitr_1.37
## [25] kableExtra_1.3.4 ggplot2_3.3.5 viridis_0.6.2
## [28] viridisLite_0.4.0 pbapply_1.5-0 readxl_1.3.1
## [31] lubridate_1.8.0 remotes_2.4.2
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 tidyselect_1.1.2 fftw_1.0-6.1
## [4] htmlwidgets_1.5.4 grid_4.1.0 munsell_0.5.0
## [7] codetools_0.2-18 DT_0.20 miniUI_0.1.1.1
## [10] withr_2.4.3 spatstat.random_2.1-0 Brobdingnag_1.2-7
## [13] colorspace_2.0-3 highr_0.9 uuid_1.1-0
## [16] rstudioapi_0.13 stats4_4.1.0 dtw_1.22-3
## [19] tensor_1.5 bayesplot_1.8.1 labeling_0.4.2
## [22] emmeans_1.8.1-1 rstan_2.21.3 polyclip_1.10-0
## [25] farver_2.1.0 bridgesampling_1.1-2 rprojroot_2.0.2
## [28] coda_0.19-4 vctrs_0.3.8 generics_0.1.2
## [31] TH.data_1.1-0 xfun_0.29 R6_2.5.1
## [34] markdown_1.1 gamm4_0.2-6 projpred_2.0.2
## [37] bitops_1.0-7 spatstat.utils_2.3-0 reshape_0.8.8
## [40] assertthat_0.2.1 promises_1.2.0.1 scales_1.1.1
## [43] multcomp_1.4-17 gtable_0.3.0 goftest_1.2-3
## [46] processx_3.5.2 xaringanExtra_0.7.0 sandwich_3.0-1
## [49] rlang_1.0.1 systemfonts_1.0.4 splines_4.1.0
## [52] spatstat.geom_2.3-2 broom_0.7.12 checkmate_2.0.0
## [55] inline_0.3.19 yaml_2.3.5 reshape2_1.4.4
## [58] abind_1.4-5 threejs_0.3.3 crosstalk_1.2.0
## [61] backports_1.4.1 httpuv_1.6.5 rsconnect_0.8.25
## [64] tensorA_0.36.2 tools_4.1.0 ellipsis_0.3.2
## [67] spatstat.core_2.4-0 raster_3.5-15 jquerylib_0.1.4
## [70] RColorBrewer_1.1-2 proxy_0.4-26 plyr_1.8.6
## [73] base64enc_0.1-3 RCurl_1.98-1.6 ps_1.6.0
## [76] prettyunits_1.1.1 rpart_4.1-15 deldir_1.0-6
## [79] zoo_1.8-9 cluster_2.1.2 magrittr_2.0.2
## [82] ggdist_3.1.0 colourpicker_1.1.1 mvtnorm_1.1-3
## [85] matrixStats_0.61.0 shinyjs_2.1.0 mime_0.12
## [88] evaluate_0.15 arrayhelpers_1.1-0 xtable_1.8-4
## [91] shinystan_2.5.0 gridExtra_2.3 rstantools_2.1.1
## [94] compiler_4.1.0 tibble_3.1.6 crayon_1.5.0
## [97] minqa_1.2.4 StanHeaders_2.21.0-7 htmltools_0.5.2
## [100] mgcv_1.8-36 later_1.3.0 RcppParallel_5.1.5
## [103] DBI_1.1.1 boot_1.3-28 permute_0.9-7
## [106] cli_3.2.0 parallel_4.1.0 igraph_1.2.11
## [109] pkgconfig_2.0.3 signal_0.7-7 terra_1.5-21
## [112] spatstat.sparse_2.1-0 xml2_1.3.3 svUnit_1.0.6
## [115] dygraphs_1.1.1.6 svglite_2.1.0 bslib_0.3.1
## [118] webshot_0.5.2 estimability_1.4.1 rvest_1.0.2
## [121] stringr_1.4.0 distributional_0.3.0 callr_3.7.0
## [124] digest_0.6.29 vegan_2.5-7 spatstat.data_2.1-2
## [127] rmarkdown_2.11 cellranger_1.1.0 shiny_1.7.1
## [130] gtools_3.9.2 rjson_0.2.21 nloptr_2.0.0
## [133] lifecycle_1.0.1 nlme_3.1-152 jsonlite_1.8.0
## [136] fansi_1.0.2 pillar_1.7.0 lattice_0.20-44
## [139] loo_2.4.1 fastmap_1.1.0 httr_1.4.2
## [142] pkgbuild_1.3.1 survival_3.2-11 glue_1.6.1
## [145] xts_0.12.1 shinythemes_1.2.0 stringi_1.7.6
## [148] sass_0.4.0